plus sn10x, SI
Baf53cre ENS neurons, SI data
CR infection 7d, CTL x4 CKO(Ahr-cko) x4
loading 60k cells,
cellranger called 35,059 cells
pure neurons downstream
load dependancies
load obj
GEX.seur <- readRDS("./GEX0112.seur.pre.rds")
GEX.seur## An object of class Seurat
## 24081 features across 18878 samples within 1 assay
## Active assay: RNA (24081 features, 1084 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapcolor.FB <- ggsci::pal_d3("category20c")(20)[c(1,6,11,16,
2,7,12,17
)]
color.cnt <- color.FB[c(1,5)]DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
filtering
GEX.seur <- subset(GEX.seur, subset = preAnno1 %in% levels(GEX.seur$preAnno1)[1:4] & DoubletFinder0.1 == "Singlet")
GEX.seur## An object of class Seurat
## 24081 features across 16562 samples within 1 assay
## Active assay: RNA (24081 features, 1084 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Tac1",
"Kctd8","Ntrk2","Penk",
"Fut9","Nfatc1","Egfr","Mgll",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))
pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s
pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s
pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s
re-clustering
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 200)## [1] "Mgat4c" "Zfp804a" "Cntn5" "Gal"
## [5] "Adgrg6" "Gm30613" "Kcnb2" "Klhl1"
## [9] "Gm32647" "Robo2" "Cntnap2" "Cdh8"
## [13] "Sst" "Cpne4" "Ntng1" "Nrxn3"
## [17] "Gm29521" "Gm21847" "Ano2" "Pdzrn4"
## [21] "Brinp3" "Lingo2" "Cdh9" "Ntrk3"
## [25] "Gpr149" "Ebf1" "Nmu" "Kcnip4"
## [29] "Nrg1" "Sgcz" "Gm15261" "Tmeff2"
## [33] "Astn2" "Dgkg" "Zfp804b" "Pcdh9"
## [37] "Cdh18" "Cmah" "Prkg1" "Clstn2"
## [41] "Car10" "Schip1" "March1" "Cdh6"
## [45] "Grm5" "Gpc5" "Nxph2" "Vip"
## [49] "Piezo2" "Chsy3" "Sema5a" "Fgf14"
## [53] "Kcnq5" "Asic2" "Vwc2l" "Rasgef1b"
## [57] "Plxna4" "Egfem1" "Dcc" "Grin3a"
## [61] "4930509J09Rik" "Unc5d" "Pbx3" "Ccbe1"
## [65] "Efr3a" "Pcdh10" "4930432L08Rik" "Itgb6"
## [69] "Cysltr2" "Zfhx3" "Spock3" "Cpa6"
## [73] "Gm20754" "Gm15680" "9530059O14Rik" "Col25a1"
## [77] "Csmd3" "Myl1" "Penk" "Pde1a"
## [81] "Kcnh7" "Robo1" "Dgkb" "Tnr"
## [85] "Serpini1" "Trhde" "Cacna2d3" "Gm26612"
## [89] "Dpp10" "Acp7" "Skap1" "Wdr17"
## [93] "Fut9" "Plcl1" "Luzp2" "Gm15584"
## [97] "Csmd1" "Abca9" "4930422I22Rik" "Nxph1"
## [101] "A330008L17Rik" "Gm49938" "9530026P05Rik" "1700034P13Rik"
## [105] "Pde4d" "Tafa2" "Gm11342" "Gm20713"
## [109] "Bmpr1b" "Tac1" "Cntn4" "Il1rapl1"
## [113] "Trps1" "Gna14" "Pcdh11x" "5530401A14Rik"
## [117] "Cux2" "Bmp5" "Cdh19" "Arhgap6"
## [121] "Gm29771" "Kcnd2" "8030451A03Rik" "Galnt18"
## [125] "Chrm3" "Kif26b" "C730002L08Rik" "Dgki"
## [129] "Gm16685" "Kctd16" "Prune2" "D030068K23Rik"
## [133] "Kctd8" "Hs6st3" "B230323A14Rik" "Rab27b"
## [137] "Myh3" "mt-Co3" "Calcb" "Ghr"
## [141] "Slc4a4" "1700111E14Rik" "Lama2" "Grik1"
## [145] "mt-Atp6" "Gm49127" "B230312C02Rik" "1700063D05Rik"
## [149] "Mme" "Iqgap2" "Ano5" "Galntl6"
## [153] "Nos1" "Gm11339" "Ttc29" "4930587E11Rik"
## [157] "Chrna9" "Olfr1369-ps1" "Gm12239" "Gm20560"
## [161] "Gm26632" "mt-Co1" "Gm16226" "Galr1"
## [165] "AC150683.1" "Gm30382" "Rasgrf1" "Gm28494"
## [169] "Fbxl7" "Gm49678" "Gm40841" "Rxfp1"
## [173] "C79798" "Rerg" "Meis2" "mt-Co2"
## [177] "Gm32916" "Grm7" "1700003I16Rik" "Adam2"
## [181] "Cadm2" "Adamts9" "Synpr" "Flrt2"
## [185] "Alcam" "Thsd7b" "Cntn3" "Nell1"
## [189] "Met" "Thsd4" "Plxdc2" "Tbx20"
## [193] "Sorcs3" "4931415C17Rik" "Dock10" "Cdh10"
## [197] "Ephb1" "Dapk2" "4833423E24Rik" "Mid1"GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))## Centering and scaling data matrixPCA
# exclude MT genes and more
#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)
DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))## PC_ 1
## Positive: Ntrk3, Ano2, Tmeff2, Robo2, Cdh8, Cpne4, Nrxn3, Mgat4c, Adgrg6, Clstn2
## Cntn5, Pdzrn4, Myl1, Gpr149, Spock3, Zfp804a, Cdh6, Cysltr2, Pcdh10, Cacna2d3
## Grin3a, Sgcz, Dgkg, Astn2, Iqgap2, Kif26b, Plxna4, Arhgap6, Cux2, Itgb6
## Negative: Lrrtm4, Galntl6, Tox, Ndst4, Grik1, Bnc2, Pcdh7, Cacnb2, Tshz2, Cdc14a
## Kcnd2, Epha5, Chrna7, Specc1, Lrp1b, Galnt17, Syt6, Plcxd3, Zfp536, Fbxw15
## Kcnab1, Synpo2, Brinp2, Nos1, St6galnac3, Lrrc4c, Rspo2, Fbxw24, Ptprt, Man2a1
## PC_ 2
## Positive: Rbfox1, Bnc2, Ptprt, Tafa1, Gpc6, Tshz2, Grik1, Frmd4b, St6galnac3, Fbxw15
## Brinp2, Tcf7l2, Cdc14a, Plcxd3, Tox, Fbxw24, Tmem132c, Dock2, Pcdh7, Pld5
## Caln1, Oprk1, Unc5c, Specc1, Adamtsl1, Slc35f4, Stxbp5l, Nfia, Ccbe1, Sdk2
## Negative: Nos1, Egfem1, Auts2, Kcnq5, Cadps2, Etv1, Gfra1, Epha5, Dach1, Asic2
## Kcnab1, Alcam, Schip1, Dgkb, Cmah, Kcnt2, Stxbp6, Rgs6, Hs6st3, Tmem108
## Adgrl3, Lrrc4c, Creb5, Cdh11, Fat3, Il1rapl1, Kcnj3, Ebf1, Tenm3, Cacnb2
## PC_ 3
## Positive: Sgcd, Adgrg6, Cysltr2, Gpr149, Nos1, Nmu, Slc4a4, Grin3a, Nfia, Dapk2
## Itgb6, Fgf13, Ano2, Ccbe1, Rgs6, Dgkg, Cbln2, Kcnab1, Lhfpl2, Ngfr
## Cpne4, Otof, Zfp536, Syt2, Gfra1, Efr3a, Zfp804a, Pex7, Calcb, Smad6
## Negative: Cdh18, Kcnip4, Csmd3, Klhl1, Kctd8, Cadm2, Pbx3, Gpc6, Htr4, Cntn3
## Meis1, Dlc1, Gabrg3, Skap1, C79798, Csmd1, Tenm2, March1, Car10, Serpini1
## Pde10a, Sema5a, Vwc2l, Khdrbs2, Dmd, Pde4d, Cdh9, Sybu, Sv2c, Pakap
## PC_ 4
## Positive: March1, Klhl1, Sema5a, Vwc2l, Cdh9, Pbx3, Kcnh7, Prune2, Zfhx3, Chsy3
## C79798, Ntng1, Rasgef1b, Tenm4, Galnt18, Bcl2, Dcc, Galnt17, Zbbx, Ebf1
## Il1rapl1, Lncbate1, Olfr78, Kcnb2, Alk, Sez6l, Scgn, Bmpr1b, Kif26b, Etv1
## Negative: Dock10, Lingo2, Kcnip4, Ndst4, Nxph1, Tac1, Gda, Sorcs1, Prkg1, Epha5
## Csmd3, Cntn5, Thsd4, Dmd, Lrrtm4, Lrrc7, Kctd8, Kcnt2, Lrp1b, Fgf13
## Cntn3, Lama2, Ccdc60, Nyap2, Hgf, Pld5, Nhs, Ryr3, Calcrl, Grem2
## PC_ 5
## Positive: Trhde, Cntn4, Ebf1, Trps1, Nrg1, Ptprd, Cpa6, Csmd1, Zmat4, Col18a1
## Kcnd2, Rmst, Egfem1, Gal, Chsy3, Myo1e, Nkain3, Npas3, Luzp2, Nav2
## Shisa6, Fstl4, Moxd1, Sctr, Prkg2, Myo16, Ntng1, Adgrl2, Asic2, Prkd1
## Negative: Dgkb, Kcnt2, Thsd7b, Klhl1, Alcam, Prkg1, Vwc2l, Rasgef1b, Il1rapl1, Fgf13
## Pbx3, Lrrc4c, Cdh20, Cdh9, Dpp10, Vcan, Khdrbs2, Zfhx3, Mgat4c, Galnt18
## Stard13, C79798, P3h2, Auts2, Zbbx, Slc44a5, Olfr78, Nmur2, Slc2a13, Rgs6length(setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene)))## [1] 1046setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene))[1:300]## [1] "Mgat4c" "Zfp804a" "Cntn5" "Gal" "Adgrg6" "Kcnb2"
## [7] "Klhl1" "Robo2" "Cntnap2" "Cdh8" "Sst" "Cpne4"
## [13] "Ntng1" "Nrxn3" "Ano2" "Pdzrn4" "Brinp3" "Lingo2"
## [19] "Cdh9" "Ntrk3" "Gpr149" "Ebf1" "Nmu" "Kcnip4"
## [25] "Nrg1" "Sgcz" "Tmeff2" "Astn2" "Dgkg" "Zfp804b"
## [31] "Pcdh9" "Cdh18" "Cmah" "Prkg1" "Clstn2" "Car10"
## [37] "Schip1" "March1" "Cdh6" "Grm5" "Gpc5" "Nxph2"
## [43] "Vip" "Piezo2" "Chsy3" "Sema5a" "Fgf14" "Kcnq5"
## [49] "Asic2" "Vwc2l" "Rasgef1b" "Plxna4" "Egfem1" "Dcc"
## [55] "Grin3a" "Unc5d" "Pbx3" "Ccbe1" "Efr3a" "Pcdh10"
## [61] "Itgb6" "Cysltr2" "Zfhx3" "Spock3" "Cpa6" "Col25a1"
## [67] "Csmd3" "Myl1" "Penk" "Pde1a" "Kcnh7" "Robo1"
## [73] "Dgkb" "Tnr" "Serpini1" "Trhde" "Cacna2d3" "Dpp10"
## [79] "Acp7" "Skap1" "Wdr17" "Fut9" "Plcl1" "Luzp2"
## [85] "Csmd1" "Abca9" "Nxph1" "Pde4d" "Tafa2" "Bmpr1b"
## [91] "Tac1" "Cntn4" "Il1rapl1" "Trps1" "Gna14" "Pcdh11x"
## [97] "Cux2" "Bmp5" "Cdh19" "Arhgap6" "Kcnd2" "Galnt18"
## [103] "Chrm3" "Kif26b" "Dgki" "Kctd16" "Prune2" "Kctd8"
## [109] "Hs6st3" "Rab27b" "Myh3" "Calcb" "Ghr" "Slc4a4"
## [115] "Lama2" "Grik1" "Mme" "Iqgap2" "Ano5" "Galntl6"
## [121] "Nos1" "Ttc29" "Chrna9" "Galr1" "Rasgrf1" "Fbxl7"
## [127] "Rxfp1" "C79798" "Rerg" "Meis2" "Grm7" "Adam2"
## [133] "Cadm2" "Adamts9" "Synpr" "Flrt2" "Alcam" "Thsd7b"
## [139] "Cntn3" "Nell1" "Met" "Thsd4" "Plxdc2" "Tbx20"
## [145] "Sorcs3" "Dock10" "Cdh10" "Ephb1" "Dapk2" "Mid1"
## [151] "Gabrg3" "Cbln2" "Tenm2" "Hs3st4" "Npas3" "Galnt13"
## [157] "Ppp3ca" "Stk31" "Mast4" "Loxl1" "Arhgap15" "Cadps2"
## [163] "Olfr78" "Ror1" "Gpm6a" "Prr16" "Dach1" "Rbpms"
## [169] "Tenm4" "Rmst" "Tenm1" "Zfp286os" "Zbbx" "Hccs"
## [175] "Ngfr" "Aff2" "Nox3" "Ntm" "Calcrl" "Wbp2nl"
## [181] "Etv1" "Sctr" "Mecom" "Hcn1" "Col8a1" "Pla2r1"
## [187] "Tafa1" "Rarb" "Otof" "Adamts6" "Nell1os" "Rtl4"
## [193] "Sorcs1" "Atg4a" "Cpne8" "Lrrc4c" "Dab1" "Vcan"
## [199] "Wee2" "Kcna1" "Nwd1" "Auts2" "Serpine2" "F13a1"
## [205] "Stxbp6" "Kirrel3" "Bfsp2" "Necab1" "Abca8a" "Itgbl1"
## [211] "Fstl4" "Lhfpl2" "Nek1" "Ptprt" "Pantr1" "Il18r1"
## [217] "Btk" "Grp" "Slc44a5" "Slc2a13" "Serpinb7" "Serpinb5"
## [223] "Nkain3" "Gabrb1" "Xlr4a" "Ptk7" "Opcml" "Rbfox1"
## [229] "Moxd1" "L3mbtl4" "Igsf21" "Rgs6" "Prkag2" "Fpr1"
## [235] "Col18a1" "Gfra1" "Lmcd1" "St18" "Brinp2" "Slc44a4"
## [241] "Edil3" "Naaladl2" "Slco1a6" "Ropn1" "Slc6a7" "Bcl2"
## [247] "Scgn" "Casp8" "Ptprd" "Grid1" "Lrp1b" "Htr4"
## [253] "Cntnap5b" "Nfatc1" "Tex15" "Aff3" "Bnc2" "Scn7a"
## [259] "Shisa9" "Gda" "Tacr1" "Inhbc" "Prkca" "Gpc6"
## [265] "Esr1" "Pakap" "Runx1" "Kcnn3" "Hs3st2" "Ntrk2"
## [271] "Apol7e" "Slco1a1" "Zbtb7c" "Siah3" "Mgat4e" "Hsd17b3"
## [277] "Lncbate1" "Eya1" "Bche" "Spag16" "Plcb1" "Adgrl2"
## [283] "Lrrc43" "Vmn2r65" "Tmprss4" "Antxr2" "Pde10a" "Scnn1a"
## [289] "Scube1" "Lin9" "Pde7b" "Nsun2" "Htr3b" "Apol7b"
## [295] "Gata4" "Olfr92" "Sv2c" "Adcy2" "Cask" "Dmd"DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 16562
## Number of edges: 641121
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8276
## Number of communities: 26
## Elapsed time: 2 secondsRun UMAP/tSNE
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 145)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 00:10:26 UMAP embedding parameters a = 0.9922 b = 1.112## 00:10:26 Read 16562 rows and found 18 numeric columns## 00:10:26 Using Annoy for neighbor search, n_neighbors = 20## 00:10:26 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 00:10:28 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpkH5UfV\file4b6c6ed3248
## 00:10:28 Searching Annoy index using 1 thread, search_k = 2000
## 00:10:32 Annoy recall = 100%
## 00:10:32 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 00:10:33 Initializing from normalized Laplacian + noise (using irlba)
## 00:10:35 Commencing optimization for 200 epochs, with 482664 positive edges
## 00:10:52 Optimization finishedDimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
length(as.vector(unlist(markers.new.s)))## [1] 92FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))
FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))
Anno
#
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(1,2,0,13,12)] <- "EMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(8)] <- "EMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9,14)] <- "EMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(18,24)] <- "EMN4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3,6,5,4)] <- "IMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(16)] <- "IMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "IMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(23)] <- "IMN4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7)] <- "IN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(19)] <- "IN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(22)] <- "IN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15,11)] <- "IPAN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(17,20)] <- "IPAN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(25)] <- "IPAN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(21)] <- "IPAN4"
GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
levels = c(paste0("EMN",1:4),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",1:4)))
#
GEX.seur$Anno2 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(1,2,0,13,12)] <- "EMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(8)] <- "EMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(9,14)] <- "EMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(18)] <- "EMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(24)] <- "EMN5"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(3,6,5,4)] <- "IMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(16)] <- "IMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(10)] <- "IMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(23)] <- "IMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(7)] <- "IN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(19)] <- "IN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(22)] <- "IN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(15)] <- "IPAN1.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(11)] <- "IPAN1.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(17)] <- "IPAN2.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(20)] <- "IPAN2.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(25)] <- "IPAN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(21)] <- "IPAN4"
GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))))DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)
cells.nn <- colnames(GEX.seur)[GEX.seur$preAnno2 %in% c("IPAN3","IPAN4") & GEX.seur$Anno1 == "IMN4"]DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "Anno1", cells.highlight = cells.nn) +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)
#
GEX.seur$Anno3 <- as.character(GEX.seur$Anno1)
GEX.seur$Anno3[cells.nn] <- "Mixlike"
GEX.seur$Anno3 <- factor(GEX.seur$Anno3,
levels = c(levels(GEX.seur$Anno1),"Mixlike"))DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno3", label.size = 3.2) +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "FB.info", cols = color.FB)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "FB.info", cols = color.FB)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno3")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "Anno3")
length(cells.nn)## [1] 42GEX.seur <- subset(GEX.seur, subset = Anno3 != "Mixlike")
GEX.seur## An object of class Seurat
## 24081 features across 16520 samples within 1 assay
## Active assay: RNA (24081 features, 1046 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno2")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "Anno2")
markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "Anno2"
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
# test.use = "MAST",
# #test.use = "wilcox",
# logfc.threshold = 0.25)
GEX.markers.pre <- read.table("Anno2.20240125.markers.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
levels = levels(GEX.seur$Anno2))
markers.pre_t32 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.1) %>%
top_n(n = 32, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t48 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.pre$gene,invert = T,value = T)) %>%
top_n(n = 48, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
filter(p_val_adj < 0.01) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
filter(p_val_adj < 0.01) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$genettt = 557
ttt/60## [1] 9.283333ttt/64## [1] 8.703125ttt/65## [1] 8.569231pp.t60 <- list()
for(i in 1:9){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(64*i-63):(64*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t60## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
ttt = 972
ttt/60## [1] 16.2ttt/64## [1] 15.1875ttt/65## [1] 14.95385pp.t120 <- list()
for(i in 1:15){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(65*i-64):(65*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
#pp.t120table(GEX.seur$DoubletFinder0.05)##
## Singlet
## 16520table(GEX.seur$DoubletFinder0.1)##
## Singlet
## 16520markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Itga6","Cntnap5b",
"Pgm5","Chgb",
"Nxph1",
"Gabrb1","Glp2r","Nebl",
"Lrrc7",
"Ryr3","Eda",
"Hgf","Lama2","Efnb2",
"Tac1",
"Kctd8",
"Ptn",
"Ntrk2","Penk","Itgb8",
"Fut9","Nfatc1","Egfr","Pdia5",
"Ahr","Mgll",
"Aff3",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2",
"Col25a1",
"Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Nova1",
"Cdh10","Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Itgb6","Met","Itgbl1",
"Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9",
"Car10","Dcc",
"Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))
pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s
pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s
pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s
cell composition
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAGTCTA-1 CR7d_SI 2107 1263 0.00000000 0.4271476
## AAACCCAAGACATACA-1 CR7d_SI 2345 1327 0.04264392 0.3411514
## AAACCCAAGACTCATC-1 CR7d_SI 5340 2565 0.71161049 0.6179775
## AAACCCAAGAGGTTTA-1 CR7d_SI 2542 1295 0.27537372 0.4327301
## AAACCCAAGGCTCTAT-1 CR7d_SI 2705 1514 1.25693161 0.7024030
## AAACCCAAGGTCTTTG-1 CR7d_SI 1053 765 0.00000000 0.5698006
## FB.info S.Score G2M.Score Phase RNA_snn_res.2
## AAACCCAAGAAGTCTA-1 CTL.3 0.007075537 -0.006547972 S 6
## AAACCCAAGACATACA-1 CTL.3 -0.002734042 0.015554703 G2M 19
## AAACCCAAGACTCATC-1 CTL.4 0.082419696 0.005693510 S 15
## AAACCCAAGAGGTTTA-1 CTL.1 -0.003888110 0.010933631 G2M 11
## AAACCCAAGGCTCTAT-1 CKO.4 0.008380733 0.021945733 G2M 11
## AAACCCAAGGTCTTTG-1 CKO.3 0.008806639 -0.010166359 S 20
## seurat_clusters cnt pANN_0.25_0.005_944 DoubletFinder0.05
## AAACCCAAGAAGTCTA-1 6 CTL 0.0000000 Singlet
## AAACCCAAGACATACA-1 19 CTL 0.1111111 Singlet
## AAACCCAAGACTCATC-1 15 CTL 0.1190476 Singlet
## AAACCCAAGAGGTTTA-1 11 CTL 0.0000000 Singlet
## AAACCCAAGGCTCTAT-1 11 CKO 0.0000000 Singlet
## AAACCCAAGGTCTTTG-1 20 CKO 0.0000000 Singlet
## pANN_0.25_0.005_1888 DoubletFinder0.1 sort_clusters preAnno1
## AAACCCAAGAAGTCTA-1 0.0000000 Singlet 7 IMN
## AAACCCAAGACATACA-1 0.0952381 Singlet 20 IN
## AAACCCAAGACTCATC-1 0.1190476 Singlet 15 IPAN
## AAACCCAAGAGGTTTA-1 0.0000000 Singlet 11 IPAN
## AAACCCAAGGCTCTAT-1 0.0000000 Singlet 11 IPAN
## AAACCCAAGGTCTTTG-1 0.0000000 Singlet 21 IPAN
## preAnno2 Anno1 Anno2 Anno3
## AAACCCAAGAAGTCTA-1 IMN1 IMN1 IMN1 IMN1
## AAACCCAAGACATACA-1 IN2 IN2 IN2 IN2
## AAACCCAAGACTCATC-1 IPAN1 IPAN1 IPAN1.1 IPAN1
## AAACCCAAGAGGTTTA-1 IPAN1 IPAN1 IPAN1.2 IPAN1
## AAACCCAAGGCTCTAT-1 IPAN1 IPAN1 IPAN1.2 IPAN1
## AAACCCAAGGTCTTTG-1 IPAN2 IPAN2 IPAN2.2 IPAN2GEX.seur$FB.info <- factor(as.character(GEX.seur$FB.info),
levels = c(paste0("CTL.",1:4),
paste0("CKO.",c(1,3,4))))GEX.seur$rep <- paste0("rep",
gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info",
ncol = 4, cols = color.FB)
#DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info",
# ncol = 3, cols = color.FB)DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
heatmap.1
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno1),
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(4,8,11),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno1)),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(4,8,11),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
barplot.1
stat_Anno1 <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]
stat_Anno1.s <- stat_Anno1 %>%
group_by(cnt,rep,Anno1) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom.Anno1 <- stat_Anno1.s %>%
ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title="Cell Composition of SI data, Anno1", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno1, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom.Anno1
sig.test
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.Anno1 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat_Anno1.s_N <- stat_Anno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat_Anno1.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat_Anno1.s_N$total <- total.N[paste0(stat_Anno1.s_N$cnt,
"_",
stat_Anno1.s_N$rep),"total"]
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat_Anno1.s_N$Anno1)){
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.Anno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno1)))
rownames(glm.nb_anova.Anno1_df) <- names(glm.nb_anova.Anno1)
colnames(glm.nb_anova.Anno1_df) <- gsub("X","C",colnames(glm.nb_anova.Anno1_df))
glm.nb_anova.Anno1_df## EMN1 EMN2 EMN3 EMN4 IMN1 IMN2 IMN3
## CTLvsCKO 0.3344298 0.5904035 0.5503359 0.4107641 0.3560548 0.1619612 0.7252822
## IMN4 IN1 IN2 IN3 IPAN1 IPAN2
## CTLvsCKO 0.7400615 0.0004540387 0.487195 0.3854356 0.08110179 0.9761208
## IPAN3 IPAN4
## CTLvsCKO 0.1452842 0.415677round(glm.nb_anova.Anno1_df,5)## EMN1 EMN2 EMN3 EMN4 IMN1 IMN2 IMN3 IMN4 IN1
## CTLvsCKO 0.33443 0.5904 0.55034 0.41076 0.35605 0.16196 0.72528 0.74006 0.00045
## IN2 IN3 IPAN1 IPAN2 IPAN3 IPAN4
## CTLvsCKO 0.4872 0.38544 0.0811 0.97612 0.14528 0.41568heatmap.2
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno2),
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno2)),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
barplot.2
stat_Anno2 <- GEX.seur@meta.data[,c("cnt","rep","Anno2")]
stat_Anno2.s <- stat_Anno2 %>%
group_by(cnt,rep,Anno2) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom.Anno2 <- stat_Anno2.s %>%
ggplot(aes(x = Anno2, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title="Cell Composition of SI data, Anno2", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno2, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom.Anno2
sig.test
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.Anno2 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat_Anno2.s_N <- stat_Anno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat_Anno2.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat_Anno2.s_N$total <- total.N[paste0(stat_Anno2.s_N$cnt,
"_",
stat_Anno2.s_N$rep),"total"]
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat_Anno2.s_N$Anno2)){
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.Anno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno2)))
rownames(glm.nb_anova.Anno2_df) <- names(glm.nb_anova.Anno2)
colnames(glm.nb_anova.Anno2_df) <- gsub("X","C",colnames(glm.nb_anova.Anno2_df))
glm.nb_anova.Anno2_df## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2
## CTLvsCKO 0.3344298 0.5904035 0.5503359 0.919722 0.05319255 0.3560548 0.1619612
## IMN3 IMN4 IN1 IN2 IN3 IPAN1.1
## CTLvsCKO 0.7252822 0.7400615 0.0004540387 0.487195 0.3854356 0.5741269
## IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## CTLvsCKO 0.01974064 0.5762993 0.4951306 0.1452842 0.415677round(glm.nb_anova.Anno2_df,5)## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4
## CTLvsCKO 0.33443 0.5904 0.55034 0.91972 0.05319 0.35605 0.16196 0.72528 0.74006
## IN1 IN2 IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## CTLvsCKO 0.00045 0.4872 0.38544 0.57413 0.01974 0.5763 0.49513 0.14528 0.41568